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ABSTRACT 

We investigate the formation of metal-free, Population III (Pop III), stars within a 
minihalo at z ^ 20 with a smoothed particle hydrodynamics (SPH) simulation, start- 
ing from cosmological initial conditions. Employing a hierarchical, zoom-in procedure, 
we achieve sufficient numerical resolution to follow the collapsing gas in the center of 
the minihalo up to number densities of 10^^ cm^''. This allows us to study the pro- 
tostellar accretion onto the initial hydrostatic core, which we represent as a growing 
sink particle, in improved physical detail. The accretion process, and in particular its 
termination, governs the final masses that were reached by the first stars. The primor- 
dial initial mass function (IMF), in turn, played an important role in determining to 
what extent the first stars drove early cosmic evolution. We continue our simulation 
for 5000 yr after the first sink particle has formed. During this time period, a disk- 
like configuration is assembled around the first protostar. The disk is gravitationally 
unstable, develops a pronounced spiral structure, and fragments into several other 
protostellar seeds. At the end of the simulation, a small multiple system has formed, 
dominated by a binary with masses ^ 40 Mq and ~ 10 Mq. If Pop III stars were to 
form typically in binaries or small multiples, the standard model of primordial star 
formation, where single, isolated stars are predicted to form in minihaloes, would have 
to be modified. This would have crucial consequences for the observational signature 
of the first stars, such as their nucleosynthetic pattern, and the gravitational-wave 
emission from possible Pop III black-hole binaries. 
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1 INTRODUCTION 

How did the cosmic dark ages end, and how was the homo- 
geneous early Universe transformed into the highly complex 
state that we observe today (e.g. Barkana & Loeb 2001; 
Miralda-Escude 2003; Ciardi & Ferrara 2005)? One of the 
key questions is to understand the formation and properties 
of the first stars, the so-called Population III (or Pop III), 
since they likely played a crucial role in driving early cosmic 
evolution (e.g. Bromm & Larson 2004; Glover 2005; Bromm 
et al. 2009). For instance, the radiation from the first stars 
contributed to the reionization of the intergalactic medium 
(IGM), at least during its initial phase (e.g. Kitayama et al. 
2004; Sokasian et al. 2004; Whalen et al. 2004; Alvarez et al. 
2006; Johnson et al. 2007). After the first stars exploded as 
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supernovae (SNe) they spread heavy elements through the 
IGM, thereby providing the initial cosmic metal enrichment 
(e.g. Madau et al. 2001; Mori et al. 2002; Bromm et al. 2003; 
Wada & Venkatesan 2003; Norman et al. 2004; Greif et al. 
2007; Tornatore et al. 2007). The resulting metallicity was 
hkely crucial in governing the transition from early high- 
mass star formation to the normal-mass star formation seen 
today (e.g. Omukai 2000; Bromm et al. 2001; Schneider et al. 
2002; Bromm & Loeb 2003; Frebel et al. 2007, 2009; Smith 
& Sigurdsson 2007). 

The extent to which the radiation and metal enrichment 
from the first stars infiuenced the early Universe is largely 
determined by their mass. Currently, Pop III stars are be- 
fieved to have first formed around z ~ 20 within minihaloes 
of mass M ~ lO'* Mq (e.g. Haiman et al. 1996; Tegmark 
et al. 1997; Yoshida et al. 2003). The potential wells of 
these minihaloes are provided by dark matter (DM), which 
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serve to initially gather primordial gas into them. This is in 
contrast to present-day star-forming giant molecular clouds 
(GMCs), which are of similar mass to minihaloes, but which 
form within much larger galaxies and are not DM dom- 
inated. Previous three-dimensional simulations (e.g. Abel 
et al. 2002; Bromm et al. 2002) have studied the forma- 
tion of M ~ 1000 Mq prestellar cores with sizes < 0.5 pc 
within such minihaloes. The density in these cores is similar 
to those in present-day ones, but having considerably higher 
temperature (~ 200 K for primordial cores versus 10 — 15 K 
for those within GMCs). 

The mass of the first stars depends on the further evo- 
lution of such cores. Similar to present-day star formation 
(e.g. McKee & Tan 2002), it is expected that some of the 
gas within these cores will gravitatioually collapse into one 
or more protostars of initial mass ^ 5 x 10~^ Mq (Omukai 
& Nishi 1998; Yoshida et al. 2008). A fraction of the core 
mass will be accreted onto the protostar(s), ultimately form- 
ing stars much more massive than the initial seeds. Three- 
dimensional simulations which studied the growth of a pri- 
mordial protostar through accretion imply that the first 
stars were quite meissive, finding an upper limit of ~ 500 Mq 
(Bromm & Loeb 2004). Some of the earliest studies of Pop 
III star formation, however, predicted that Pop III stellar 
masses might be quite low, < 1 Mq (e.g. Kashlinsky & Rees 
1983). These authors emphasized the importance of angu- 
lar momentum in determining the mass of Pop III stars, 
predicting that rotational effects would cause the primor- 
dial gas clouds to collapse into a dense disk. Only after the 
disk cooled to ~ 1000 K through H2 line emission could 
fragmentation occur. Depending on the spin of the cloud, 
the fragments thus formed could have had masses as low as 
0.2 M0. Later studies by Saigo et al. (2004) also found 
that high initial angular momentum would allow primordial 
clouds to form disks, and they further predicted that disk 
fragmention into binaries would bo common. Similar calcu- 
lations by Machida et al. (2008), which wore initialized at 
much lower densities and therefore at an earlier point in 
the gas collapse, found fragmentation would occur even for 
primordial gas clouds with very modest initial angular mo- 
mentum. Finally, the recent work by Clark et al. (2008) in- 
cluded a simulation of a 500 Mq primordial gas cloud which 
fragmented into ~ 20 stellar objects. However, none of the 
above-mentioned investigations were initialized on cosmo- 
logical scales. 

Our goal in this study is to further improve our un- 
derstanding of how a Pop III star is assembled through the 
process of accretion within the first minihaloes. More specif- 
ically, we aim to determine the final masses reached by the 
first stars, and whether the gas within a host minihalo forms 
a single star or fragments into a binary or small multiple 
system. To this end we perform a three-dimensional cosmo- 
logical simulation that follows the formation of a minihalo 
with very high numerical resolution, allowing us to trace 
the evolution of its gas up to densities of 10^"^ cm^'\ Un- 
like the study of Bromm & Loeb (2004), which employed 
a constrained initialization technique on scales of ~ 150 pc, 
our simulation is initialized on cosmological scales. This al- 
lows the angular momentum of the mirnhalo to be found 
self-consistently, eliminating the need for an assumed spin 
parameter. In a fashion similar to Bromm & Loeb (2004), 
we represent high-density peaks that are undergoing gravita^ 



tional runaway collapse using the sink particle method first 
introduced by Bate et al. (1995), allowing us to follow the 
mass flow onto the sinks well beyond the onset of collapse. 

The recent simulation by Turk et al. (2009), starting 
from cosmological initial conditions, has shown fragmenta- 
tion and the formation of two high-density peaks, which are 
interpreted as seeds for a Pop III binary. These peaks axe 
similar to the multiple sinks found in our simulation. Al- 
though Turk et al. employed even higher numerical reso- 
lution, they were able to follow the evolution of the high- 
density peaks for only ^ 200 yr. Our sink particle method, 
on the other hand, enabled us to continue our calculation 
for several thousand years after initial runaway collapse. We 
could thus follow the subsequent disk evolution and further 
fragmentation in detail, rendering our study ideally comple- 
mentary to their work. 

The structure of our paper is as follows. In Section 2, we 
discuss the basic physical ingredients determining the accre- 
tion flow, whereas we describe our numerical methodology 
in Section 3. We proceed to present our results (Section 4), 
and conclude in Section 5. 



2 PHYSICAL INGREDIENTS 
2.1 Protostellar accretion 

Unlike the gas from which present-day stars are formed, 
which can typically cool to ~ 10 K, primordial gas contained 
no metals and could only cool down to ~ 100 K through H2 
and HD line transitions. These elevated temperatures lead 
to a higher accretion rate, which can be estimated on dimen- 
sional grounds by assuming that the Jeans mass is infalling 
at the free-fall rate, resulting in 
3 

M ~ ^ oc T^/^ , (1) 

Cr 

where Ca is the soundspeed. This is quite close to the Shu 
(1977) similarity solution for collapse of a singular isother- 
mal sphere, which gives M = 0.9754 /G. The Shu solution 
begins with an isothermal cloud with p cx r^^ and irutially 
negligible infall velocities. The solution then describes the 
subsequent 'inside-out' collapse of the cloud, which physi- 
cally corresponds to accretion of mass onto the central pro- 
tostar. The earlier similarity solution of Larson (1969) and 
Pcuston (1969) described the collapse of an initially urn- 
form cloud into a peaked structure with an density pro- 
file. This solution was later expanded by Hunter (1977) to 
capture the evolution after the initial formation of a hy- 
drostatic core when it grows via accretion of the envelope. 
At the point of initial protostar formation, corresponding to 
when the central density becomes infinite in the similarity 
solution, the infall velocity is already 3.3 times the sound 
speed, and the accretion rate grows to 48 times higher than 
that found in Shu collapse (Hunter 1977). 

In our simulation, the collapse of the gas from initially 
more uniform densities can be expected to behave similarly 
to a Larson-Penston-type solution (e.g. Omukai & Nishi 
1998). The following accretion of this gas onto the cen- 
tral hydrostatic core should proceed from the 'inside-out' at 
rates similar to those found in the Shu and Larson- Penston- 
Hunter solutions. 
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Figure 1. Final simulation output shown in density projection along the x-z plane on progressively smaller scales, as labeled on each 
panel. White boxes denote the region to be depicted in the next-smallest scale. In the bottom two panels, the asterisk denotes the 
location of the first sink formed. Top left: Entire simulation box. Top right: Minihalo and surrounding filamentary structure. Bottom 
right: Central fOO pc of minihalo. Bottom left: Central 10 pc of minihalo. Note how the morphology approaches an increasingly smooth, 
roughly spherical distribution on the smallest scales, where the gas is in quasi-hydrostatic equilibrium, just prior to the onset of runaway 
collapse. 



2.2 Protostellar feedback effects 



As the protostellar mass grows through accretion, feedback 
from photoionization, heating, and radiation pressure will 
begin to have important effects on the accretion flow and 
may set the upper mass limit for both present-day and pri- 
mordial stars. One important difference between feedback in 
primordial and present-day massive star formation is due to 
the existence of dust grains in the latter case. For instance, 
Ripamonti et al. (2002) found with their one-dimensional 
calculations that Pop III protostellar radiation pressure will 
not affect the initial collapse and early stages of accretion 
due to the low opacity of the infalling gas. In contrast, for 
present-day star formation the increased opacity from dust 
enhances the radiation pressure and can halt infall onto stars 
of more than a few tens of solar masses without further 
mechanisms to alleviate the radiative force (see summary in 
McKee & Ostriker 2007). On the other hand, later in the 
life of a present-day star, dust absorption of ionizing pho- 
tons will alter the evolution of the HII region around the 
star, thus reducing the photoionization feedback. 



Omukai & Palla (2003) performed a detailed study of 
spherical accretion of metal-free gas onto a growing hy- 
drostatic core. In this case, electron scattering is the ma- 
jor source of opacity determining the Eddington luminosity 
Ledd and the critical accretion rate Afcrit, above which the 
protostellar luminosity exceeds Ledd before the protostar 
reaches the zero-age main sequence (ZAMS), disrupting fur- 
ther accretion. As estimated in Bromm & Loeb (2004), if we 
assume the protostellar luminosity before the onset of hy- 
drogen burning is approximately the accretion luminosity. 
Lace — GMfM / R-M, and set this luminosity equal to Ledd, 
we find that 

M.,.^^^^^-^>clO-'M^yr-\ (2) 

where 7?* ~ SR© (Bromm et al. 2001). This approximation 
is quite close to that found in the calculations of Omukai & 
Palla (2003). It should be noted that initial accretion rates 
higher than Merit are possible since protostars begin with 
much larger radii and only gradually shrink during Kelvin- 
Helmholtz contraction. 
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Figure 2. Physical state of the collapsing core. The situation is shown just prior to the formation of the first sink particle, (a) Temperature; 
(b) free electron fraction; (c) molecular hydrogen fraction; and (d) adiabatic exponent 7^^^ vs. density. The gas that collapses into the 
minihalo is heated adiabatically to ~ 1000 K as it reaches n ~ 1 cm~^ (panel a). At this point becomes high enough (panel c) 
to allow the gas to cool through H2 rovibrational transitions to a minimum of T ~ 200 K, which is reached at n ~ 10* cm ^ (panel 
a) . After the onset of gravitational instability, the gas reaches yet higher densities, and is heated to a maximum of ~ 1000 K once it 
approaches a density of ~ 10* cm~^. At this density three-body reactions become important, and the gas is rapidly converted into fully 
molecular form (panel c). The increased H2 cooling rate then roughly equals the adiabatic heating rate, and thus the highest density gas 
is nearly isothermal. As the gas is transformed into molecular form, the equation of state changes as well so that 7ad evolves from 5/3 
to approximately 7/5 (panel d). 



Other types of feedback also have the potential to halt 
protostellar accretion and set the upper mass limit of Pop III 
stars. McKoo & Tan (2008) consider several feedback mecha- 
nisms operating during Pop III star formation, including H2 
photodissociation and Lya radiation pressure. They find, 
however, that the accretion rate is not significantly reduced 
through feedback until HII region breakout beyond the grav- 
itational escape radius, generally occuring once the star has 
reached 50-100 M©. Omukai & Inutsuka (2002), on the other 
hand, find that for spherically symmetric accretion onto an 
ionizing star, the formation of an HII region in fact does 
not impose an upper mass limit to Pop III stars. As McKee 
& Tan (2008) determine, however, this does not remain the 
case for rotating inflow. 

This previous work shows that further considerations 
such as a time-dependent accretion rate and a non-spherical, 
disk-like geometry around the protostar (see also Tan & Mc- 



Kee 2004) will be integral to fully determining the effects of 
protostellar feedback and will require study through three- 
dimensional calculations. We do not include any radiative 
feedback efi^ects here, and thus our work represents a no- 
feedback limit of protostellar accretion onto a Pop III star. 
However, during the initial stages of protostellar growth, the 
feedback should be much weaker than in later stages when 
the star becomes more massive and luminous. As mentioned 
above, McKee & Tan (2008) estimate that the accretion rate 
is not significantly reduced until the Pop HI star has gained 
~ 50 M©, and the rate of protostellar mass growth given 
by Bromm & Loeb (2004) indicates that such masses are 
not reached until ~ 10** yr after accretion begins. We thus 
end our calculation somewhat before this, after 5000 yr of 
accretion, before feedback effects have become too strong 
to ignore. Up until this point, our calculations should still 
give a good estimate of the geometry and possible fragmen- 
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Figure 3. Initial collapse prior to the formation of the first sink particle. Shown are a number of variables as a function of radial distance 
from the density maximum, (a) Number density; (b) radial velocity; (c) molecule fraction; and (d) temperature vs. r. Also shown are the 
number density and radial velocity of a Larson-Penston type similarity solution for the equation of state P = Kp'' , where 7 = 1.09 and 
K = i.2x 10" (in cgs units), as in Omukai & Nishi (1998). The solid red lines (in panels a and b) represent the latest stage of collapse 
that was calculated, when the peak density is near 10^^ cm~^, while the dashed blue lines (in panel a) represent various earlier stages 
of collapse. Note that the analytical solution reproduces the density profile found in the numerical calculation very well. The velocity 
profile is also in reasonable agreement with the simulation, though, also as in Omukai & Nishi (1998), exhibiting a somewhat larger 
normalization due to differences in initial and boundary conditions. 



tation of the prestellar core as well as the early protostellar 
accretion rate and its time dependence. 



H + H + Ha 



H2 + H2 



(4) 



2.3 Chemistry, heating and cooHng 

The chemistry, heating and cooling of the primordial gas is 
treated in a fashion very similar to previous studies (e.g. 
Bromm & Loeb 2004; Yoshida et al. 2006). We follow the 
abundance evolution of H, H"*", H~, H2, Hj, He, He^, He'''^, 
e~, and the deuterium species D, D"'", D~, HD, and HD+. 
We use the same chemical network as used in Johnson & 
Bromm (2006) and include the same cooling terms. 

At densities greater than ~ 10* cm""^, three-body pro- 
cesses gain importance. The reaction rates for 



H + H + H^Ha + H 
and 



(3) 



must be included (Palla et al. 1983). Because the gas is be- 
coming fully molecular, in our cooling rates we furthermore 
account for enhanced cooling due to collisions between Ha 
molecules. The heating rate due to Ha formation heating is 
also included. Furthermore, the equation of state of the gas 
must be modified. More specifically, we adapt the adiabatic 
exponent 7ad and the mean molecular weight fj, to account 
for the increasing Ha abundance /h2- These modifications 
for densities greater than 10* cm~^ are applied in the same 
way as described in Bromm & Loeb (2004), except for the 
modification of 7ad, which is handled in the same way as 
decribed in Yoshida et al. (2006). 
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Figure 4. Kinematics around the main sink particle, (a) Radial velocity versus distance from the main sink at 250 yr after the initial 
sink formation. Colored particles are those with densities greater than 10* cm~^ - yellow diamonds are high-density heated particles 
with temperatures greater than 2500 K, and blue circles are high-density cool particles with temperatures less than 2500 K. The red 
line shows the gravitational free-fall velocity vg, determined using the enclosed mass within the given radius. For comparison to the 
physical situation at earlier times, the green line is the radially averaged radial velocity around the main sink immediately after it first 
forms, (b) Rotational velocity versus distance from the main sink at 250 yr after sink formation. Blue and yellow particles have the same 
meaning as in panel (a). The red line shows i)Kepi calculated using the enclosed mass. The green line is the radially averaged rotational 
velocity around the main sink immediately after it first forms, (c) Radial velocity versus distance from main sink at 5000 yr after initial 
sink formation. Notation is the same as in panel (a), (d) Rotational velocity versus distance from main sink at 5000 yr after initial sink 
formation. We adopt the same convention as in panel (b). Towards later times, the cold gas settles into a nearly Keplerian configuration. 



3 NUMERICAL METHODOLOGY 
3.1 Initial setup 

We carry out our investigation using Gadget, a widely-tested 
three dimensional smoothed particle hydrodynamics (SPH) 
code (Springel et al. 2001; Springel & Hernquist 2002). 
Simulations are performed in a periodic box with size of 
100 hr^ kpc (comoving) and initialized at 2 = 99 with both 
DM particles and SPH gas particles. This is done in accor- 
dance with a ACDM cosmology with f^A = 0.7, Q,m = 0.3, 
JIb = 0.04, and h = 0.7. We use an artificially large nor- 
malization of the power spectrum, (Tg = 1.4, to accelerate 
structure formation in our relatively small box. As discussed 
in Section 5, our high value of erg is expected to have little 
effect on our overall qualitative results, though this will be 
explored further in future work. 



A preliminary run with a modest resolution of A'^sph = 
A'dm = 128"^ allows us to locate the formation site of the first 
minihalo. We then employ a standard hierarchical zoom-in 
procedure to achieve the necessary high-mass resolution in 
this region (e.g., Navarro & White 1994; Tormen et al. 1997; 
Gao et al. 2005). This method was also used in the work of 
Greif et al. (2008), in which they simulated the formation of 
a first galaxy beginning from cosmological initial conditions. 
Applying this method for our calculation involved beginning 
the simulation again ai z — 99, but now with an additional 
three nested refinement levels of length 40, 30, and 20 kpc 
(comoving) centered on the site where the first minihalo will 
form. This high-resolution run is thus initialized at high red- 
shift with the same initial fluctuations as in the preliminary 
low-resolution run, though the high-resolution region further 
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includes fluctuations at smaller wavelength scales as allowed 
by the increased Nyquist frequency. 

Particle positions are shifted so that the high resolu- 
tion region will be in the center of the simulation box. Each 
level of higher refinement replaces particles from the lower 
level with eight child particles such that in the final simula- 
tion a parent particle is replaced by up to 512 child parti- 
cles. The highest-resolution gas particles each have a mass 
msPH = 0.015 M0, so that the mass resolution of the re- 
fined simulation is; Afros — 1.5A'"ncigh'"'isPH ^ 1 M©, where 
A'neigh — 32 is the typical number of particles in the SPH 
smoothing kernel (e.g. Bate & Burkert 1997). 

The size of the refinement levels had to be chosen care- 
fully. The level of highest refinement must be large enough 
to encompass all the mass that will eventually become the 
first minihalo. If this level of highest refinement is too small, 
lower-resolution particles will drift into the center of the col- 
lapsing region and disrupt the collapse of the high-resolution 
particles. Furthermore, there will be some artificial density 
cuhanccment for lower-resolution particles that are near the 
boundary of a higher-resolution region. This is because the 
particles of lower refinement levels have higher masses, and 
for a given density the corresponding smoothing length for 
lower-resolution particles should be larger than the smooth- 
ing length for higher-resolution particles. Smoothing lengths 
for particles of both resolutions will be very similar, how- 
ever, if they are very close to each other. This causes low- 
resolution particles near a boundary to have artificially re- 
duced smoothing lengths and overly enhanced densities. 
However, this will not pose a problem as long as the highest 
refinement level is large enough such that the low-resolution 
particles are well outside of the collapsing region of interest. 



3.2 Sink particle creation 

When the density of a gas particle becomes greater than 

Wmax ~ 10^^ cm~^, this particle and the surrounding gas 
particles within a distance of less than the accretion ra- 
dius race are converted into a single sink particle. We set 
^acc equal to the resolution length of the simulation, race = 
-'-'res ~ 50 AU, where: 

L.s^0.5f^y'\ 

with />,nnx — f maxTO-H. Tlic siuk particlc's initial mass is close 
to the simulation's resolution mass, Mres — 0.7 Mq. 

Our density criterion for sink particle formation is ro- 
bust even though we do not explicitly check that the sink 
region satisfied further criteria such as being gravitationally 
bound and not having rotational support. As discussed in 
Bromm et al. (2002), crossing the Timax threshold is already 
equivalent to meeting these further criteria due to the run- 
away nature of the collapse. The average density of the disk 
(see Sec. 4.2.1) does not go above ~ 10^" cm~^, and gas must 
collapse to densities two orders of magnitude higher before 
a sink is formed at the center of the infalling region. This 
along with our very small value for race ensures that the par- 
ticles that become sinks are indeed undergoing gravitational 
collapse. Furthermore, an examination of movies of the sim- 
ulation shows that even the lowest-mass sinks are created 
only in pre-existing high-density clumps of gas that formed 



through gravitational instability, even though the clumps 
may no longer be clearly visible after the sink forms and 
accretes mass. During earlier test runs, however, using only 
the density criterion was problematic for other numerical 
reasons. This was because the larger sinks would contribute 
their high masses to the kernels of nearby gas particles, caus- 
ing them to artificially jump to n > 10^^ cm~^ and become 
sink particles themselves. These artificial sink particles had 
masses well below Mres. However, once we corrected this 
problem by removing all sink particles from the density cal- 
culation, all sinks that were formed had initial masses > 

Mres- 

After a sink's formation, its density, temperature, and 
pressure are no longer evolved in time. Instead, their density 
is held constant at 10^^ cm~^. The temperature of the sinks 
is also held constant at 650 K, a typical temperature for gas 
reaching 10^^ cm^^, and the sinks are furthermore kept at 
the pressure which corresponds to the above temperature 
and density. Giving the sink a temperature and pressure 
makes it behave more like a gas particle, and it furthermore 
prevents the existence of a pressure vacuum around the sink 
that would make the sink accretion rate artificially high (see 
Bromm et al. 2002; Martel et al. 2006). Changes in the sink's 
position, velocity, and acceleration due to gravitational in- 
teraction are still evolved, however, and as the sink accretes 
mass its gravity becomes the dominant force, and it gradu- 
ally begins to act more like a non-gaseous N-body particle. 

Once a sink particle is formed, all gas particles or other 
sink particles that cross within race of the sink particle arc 
removed from the simulation, and their mass is added to that 
of the sink particle. After each subsequent accretion event, 
as well as when the sink is initially formed, the position and 
velocity of the sink are set to be the mass- weighted average 
of the sink particle and the gas particles that it has accreted. 
The sink approximately represents the growing hydrostatic 
core, though its size, of order race, is much larger than the 
true expected size of a core (e.g. Omukai & Nishi 1998). 

The sink particle method has multiple advantages. It 
eliminates the need to include high density physics that be- 
comes relevant only at n > 10^^ cm^"'. Furthermore, simu- 
lation timesteps become prohibitively small as the density 
increases, and representing a region of n > 10^^ cm~^ with a 
single sink particle allows the mass flow into protostcUar re- 
gions to be followed for many dynamical times without con- 
tinuing to increasingly high densities and small timesteps. 
Instead of estimating the accretion rate from the instanta- 
neous density and velocity profiles at the end of the sim- 
ulation, as done in Abel et al. (2002) and Yoshida et al. 
(2006), the sink particle method allows the accretion his- 
tory to be directly followed. Furthermore, without the sink 
particle method, finding multiple fragments would require 
that they form at very nearly the same time, such as the 
fragments in Turk et al. (2009) which formed only 60 years 
apart. Using sink particles allows for the study of fragmen- 
tation that occurs on significantly longer timescales than 
this. 
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Figure 5. Density and temperature projections of the central 5000 AU. Each row shows the projections at 1000 yr (left), 2000 yr (center), 
and 5000 yr (right) after the initial sink formation. Asterisks denote the location of the most massive sink. Crosses show the location of 
the second most massive sink. Diamonds are the locations of the other sinks. Top row : Density structure of the central region in the x-y 
plane. Second row : Density structure of the central region in the x-z plane. Third row : Temperature structure of the central region in 
the x-y plane. Bottom row : Temperature structure of the central region in the x-z plane. 
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4 RESULTS 

4.1 Initial runaway collapse 

Fig. 1 shows the final simulation output on various scales. 
The expected filamentary structure is apparent, and the first 
minihalo has formed in the central region where several fila- 
ments intersect by « ~ 20. We measure the halo virialization 
radius rvir as the point at which poM = Pvir = 200pb, where 
Pdm is the average halo DM density, pvir is the DM density 
at the point of virializatiou, and pb is the average back- 
ground matter density at the time of halo virialization. We 
then find that rvir ^ 90 pc and Mhaio ^ 2.6 x 10^ M©. 
The spin parameter A — J\E\^^'^ /GM^^'^ is approximately 
0.13 for our halo, where J is the total angular momentum, 
E is the binding energy, and M is the total halo mass. For 
our minihalo A is higher than the typical value of 0.05 (e.g. 
Barnes & Efstathiou 1987; Jang-Condell & Hernquist 2001). 
This large spin may result from our enhanced value of as. 
However, as discussed in Section 5, the angular momentum 
profile of the central gas in our simulation is still quite sim- 
ilar to that found from previous cosmologicaJ simulations 
initialized with a lower as- Thus, our overall results con- 
cerning disk evolution and fragmentation are not expected 
to significantly change with a lower value of minihalo spin 

or <T8- 

After the DM in the minihalo virializes, the gas con- 
tinues to collapse. As shown in Fig. 2, the results of our 
calculation up to the point of creation of the first sink par- 
ticle agree well with previous studies (e.g. Bromm & Larson 
2004; Yoshida et al. 2006). As expected, the gas heats adia- 
batically until reaching n ^ 1 cm"'', where it attains max- 
imum temperatures of ~ 1000 K. After this point the gas 
begins to cool through H2 rovibrational transitions, reach- 
ing a minimum temperature of T 200 K at a density of 
n ~ lO'' cm~^. This is the density at which the gas reaches 
a quasi-hydrostatic 'loitering phase' (Bromm et al. 2002). 
After this phase the gas temperature rises again due to 
compressional heating until n ~ 10® cm~^, when three- 
body processes become important and cause the gas to turn 
fully molecular by n lO'^^ cm~^ (see detailed discussion in 
Yoshida et al. 2006, 2008). 

Fig. 3 shows the radial structure of the gas surround- 
ing the density pealc just before the first sink particle is 
created (compare with figure 2 in Bromm & Loeb 2004). 
The gas inside the central 100 AU has become almost fully 
molecular, and it evolves roughly isotliermally due to the 
approximate balance between compressional and H2 forma- 
tion heating and enhanced H2 cooling. The deviation from 
spherical symmetry of this central region is apparent in the 
'thickness' of the radial profiles, particularly the profile of 
the H2 fraction fy_^. However, the gas density and velocity 
profiles still have the same general properties as those of the 
spherically symmetric Larson-Penston type similarity solu- 
tions, an example of which is shown in Fig. 3. The good 
agreement is evident, at least for the initial collapse phase. 
Omukai & Nishi (1998) used the same solution to describe 
their one-dimensional radiation-hydrodynamics simulations 
of primordial protostar formation. Though the velocities in 
the similarity solution are normalized to higher values, the 
velocites in our three-dimensional calculation exhibit very 
similar overall behavior, and the density profile matches the 
similarity solution especially well. 



4.2 Protostellar accretion 

4^.2. 1 Disk formation 

We next discuss the transition from the quasi-spherical ini- 
tial collapse towards a disk-like configuration. This transi- 
tion is clearly visible in the velocity structure of the cen- 
tral region, which changes dramatically over the last sev- 
eral thousand years of the simulation as the disk-like nature 
of the central region becomes more apparent over time. In 
Fig. 4, we illustrate the evolution of the radial (left pan- 
els) and rotational velocity (right panels), measured with 
respect to the most massive sink. In comparing to the free- 
fall velocity, we calculate va using the enclosed mass Menc 
within a radial distance r from the main sink. At both times 
shown, 250 and 5000 yr after initial sink formation, pressure 
support from the central gas leads to a deviation from pure 
free fall, but the overall behavior of the radial velocity of 
the outer particles is similar to that of the vs profile. This 
is not the case for the particles very near the sink, however. 
The average radial velocity of the central gas particles, par- 
ticularly the cold ones (shown in blue), approaches zero as 
the gas begins to exhibit more rotational motion. 

The central region is evolving into a Keplerian disk 
(vrot oc r~^^^) by just 250 yr after the first sink has formed. 
Note the comparison between the particle rotational veloci- 
ties Vrot and the Keplerian velocity WKep (red line in Fig. 4), 
defined as 



f Kep = 




From Fig. 4 (panel b), it is apparent that the inner 100- 
200 AU of the disk arc becoming rotationally supported 
(vrot ^ VKep) by 250 yr after sink formation. Despite some 
amount of rotational support, however, a fraction of the par- 
ticles is still able to approach the sink, thus allowing con- 
tinued accretion. After 5000 yr, the disk has grown to a 
radial size of slightly less than 2000 AU (panel c), while the 
cold disk particles exhibit nearly Keplerian motion (panel 
d). The spike at 700 AU (panel d) corresponds to a group 
of particles rotating around the second largest sink. 

The cissembly of the disk-like configuration around the 
first sink particle, and the subsequent fragmentation into 
additional sinks, are shown in Fig. 5. The growing disk is 
unstable and develops a pronounced spiral pattern, lead- 
ing in turn to the formation of multiple high-density peaks. 
As is evident, regions of high density correlate with those 
of low temperature, due to the enhanced cooling there. The 
expanding, hour-glass shaped, bubble of hot gas (clearly vis- 
ible in the bottom row) reflects the release of gravitational 
energy when particles fall into the deepening potential well 
around the central sink. The corresponding virial tempera- 
ture is 

rp ^ GMsinkmn 4 

Ivir — -. — iU l\, 

/vB^acc 

where Mgink ~ 10 M©. This gravitational source of heating 
causes the emergence of a hot gas phase around the central 
sink (T ~ 7, 000 K), even in the absence of any direct radia- 
tive feedback from the protostar which is not included here. 
In Fig. 6, we illustrate the corresponding velocity field in 
two orthogonal planes at 5000 yr. The small infall velocities 
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onto the disk and the high rotational, roughly Keplerian, 
velocites within the disk are apparent. 

Fig. 7 shows the growth of the disk over time, where 
t = yr marks the point at which the first sink forms. We 
estimate the disk mass by including all gas particles, sinks 
excluded, with a density greater than 10* cm~^ and a spe- 
cific angular momentum within 50% of jKop, the specific 
angular momentum expected for a centrifugally supported 
disk, where jKcp = VKcpr. We here refer all velocities to the 
center of mass of the high-density (n > 10* cm~"^) parti- 
cles. In addition, we track the evolution of the mass in the 
'cool' (dotted line) and 'hot' gas phase (dashed line), corre- 
sponding to temperatures less than 2500 K and greater than 
2500 K, respectively. The mass of heated particles grows 
soon after the first sink is in place, whereas the cold compo- 
nent decreases as its mass is incorporated into the sinks or 
becomes part of the hot phase. Eventually what remains of 
the cold component is confined to the disk itself (see bottom 
two rows of Fig. 5). 

The characteristics of the star-disk system in our simu- 
lation can be compared to the analytic predictions of Tan & 
McKee (2004). At a density and temperature of 10* cm~^ 
and 1000 K, respectively, we find a value of K' ~ 1.33 from 
their equation (7), where K' is the entropy parameter, de- 
fined by 



K' = 



P/P' 



X lO^^cgs 



V300K 



lO^cm- 



(6) 



We estimate that the total star-disk system mass (red line in 
Fig. 7) is ~ 100 Mq at the end of our simulation and set 
equal to 1 for the no feedback case, where e*d is the fraction 
of infalling mass that is able to reach the star-disk system 
without being diverted by protostellar outflows, again in the 
terminology of Tan & McKee. From their equation (8), we 
then calculate an expected accretion rate onto the star-disk 
system of 10~^ Mq yr~^. In our simulation we find that 
over the 5000 yr of sink accretion, mass flows onto the star- 
disk system very steadily at approximately the same rate 
of ~ 10~^ Mq yr"'^, which is also near the value for the 
average accretion rate onto the main sink particle (see Sec. 
4.2.4). Fig. 7 further shows that once the first sink forms, the 
disk mass stays fairly constant at 35 Mq as mass falling 
from the outer regions onto the disk replaces the material 
accreted onto the sinks at similar rates. 

Tan & McKee (2004) furthermore predict the time 
needed to build up a given stellar mass in their equation 
(9). If we assume a value of 1 for /d, the ratio of disk to stel- 
lar mass, then the time to build up a 43 Mq star (the final 
mass of the main sink - see Sec. 4.2.4) is roughly 9000 yr. 
The simulation yields a similar though smaller value of 5000 
yr, which reflects the slightly higher accretion rate of our 
simulation. Finally, if we estimate a value of 1 for /kgp, the 
ratio of rotational to Keplerian velocity at the sonic point, 
then from equation (15) of Tan & McKee (2004) we find a 
predicted disk radius of ~ 3000 AU. This value is also similar 
to the value of ~ 2000 AU found from the simulation. 



4-2.2 Angular momentum evolution 

How does the distribution of angular momentum change 
during the simulation? In Fig. 8, we demonstrate the in- 
crease in rotational support of the central region over time. 



Figure 6. Velocity field of the central gas distribution. Sinks are 
denoted as in Fig. 5. Top : Density projection in the x-z plane 
after 5000 yr, shown together with the velocity field. Velocities 
are measured with respect to the center of mass of all n > 10* 
cm~'^ particles. Bottom : Same as what is shown in the top panel 
except in the orthogonal (x-y) plane. It is evident how an ordered, 
nearly Keplerian velocity structure has been established within 
the disk. 



Shown are the evolution of the enclosed mass (top panel), 
the specific angular momentum along the z-axis, (middle 
panel), and the degree of rotational support (bottom panel). 
Here, rotational support is defined as 



/rot — 



jKep 



(7) 



We determine jz by averaging over all cool (T < 2500 K) 
particles within small radial bins covering a total distance 
of lO'^AU. Note how the inner 200 Mq , corresponding to 
the central 10* AU, have spun up significantly in the pro- 
cess of disk formation. This occurs as the inner mass shells 
fall towards the center, leading to the increase in rotational 
support. The inner 200 AU of cool gas are fully rotationally 
supported by the end of the simulation (bottom panel). 
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Figure 7. Disk mass vs. time since the creation of the first sink 
particle (soUd line). Disk particles are considered to be all gas 
particles, excluding sinks, with density greater than 10** cm^"* 
and a specific angular momentum within 50% of jKcp- Once the 
first sink forms at i = yr, the red line shows the total mass 
of the star-disk system. Also shown is the total mass of n > 10® 
cm~'^ particles (sinks excluded) with temperatures less than 2500 
K (dotted line) and temperatures greater than 2500 K (dashed 
line). 



To confirm that numerical viscosity is not a dominant 
factor in determining our results, we have compared the 
torque due to numerical viscosity (tvIsc) to those exerted 
by gravity and pressure (-fgrav and Tpres). The total torque 
on a given particle within a gas cloud is given by 



— 7~grav "t" 7~prcs -1" Tvisc 

— TTlpart^part X (flgrav -t" 'Jprcs ~\~ ^^visc) , 



(8) 



where fpart is the distance of the particle from the center of 
the cloud, mpart the particle mass, and a. 



grav , (-^prcs 



and avisc 

are the accelerations due to gravity, pressure and viscosity, 
respectively. Similar to Yoshida et al. (2006) , we store these 
accelerations for each particle and every timestep. We find 
that torques from gravity and pressure are dominant by an 
order of magnitude except in regions near the sinks. The 
viscous torques within ~ 100 AU of the sinks may exceed 
the gravitational and pressure torques at any given instant. 
However, unlike -fgrav and fprcs, ivisc near the sinks does not 
tend to continuously act in any preferred direction and thus 
averages out to be much smaller than the gravitational and 
pressure torques. 



4-2.3 Binary and multiple formation 

Around 300 yr after the first sink has formed, the disk be- 
comes unstable to fragmentation and a second sink is cre- 
ated. We can estimate the value of Q in the Toomre criterion 
for disk gravitational instability, 
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Figure 8. Angular momentum structure. Top : Enclosed mass 
versus distance from main sink particle at last simulation out- 
put before sink formation (solid line) and after 5000 yr of sink 
accretion (dashed line). Middle : Specific angular momentum j^, 
of the cool (T < 2500 K) gas with respect to distance from the 
main sink at the same corresponding times. Bottom : Degree of 
rotational support of the cool gas. The solid red line shows where 
j/jKcp = li the requirement for complete rotational support. 
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tform [yr] 


Mfinal [Mq] 


rinit [AU] 
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300 


13 


60 


700 
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3700 


1.3 


930 


1110 


4 


3750 


0.8 


740 


890 


5 


4400 


1.1 


270 


240 



Q 



< 1 



(9) 



Table 1. Formation times, final masses, and distances from the 
main sink. We include all sinks still present at the end of the 
simulation. 



Here, E is the disk surface density and k is the epicyclic fre- 
quency, which is equal to the angular velocity for Keplerian 
rotation. At this point in the simulation, the soundspeed Cs 
is approximately 3 km s~^ for a ~ 1000 K disk, and k ~ 3 x 
10~^^ s~^ for a disk with Keplerian rotation at the disk outer 
radius, which we take to be i?disk — 1000 AU. The disk mass 
is close to 30 Mq, and E ~ Mdisk/vriidisk ~ 100 g/cm^. This 
gives Q ~ 0.4, well-satisfying the fragmentation criterion. 

Soon after the disk instability sets in, a pronounced spi- 
ral structure develops, as is evident in Fig. 5. The spiral per- 
turbation serves to transport angular momentum outward 
through gravitational torques, allowing some of the inner gas 
to accrete onto the sinks (see Krumholz et al. 2007). By 1000 
yr after initial sink formation, several new fragments have 
formed, represented by new sink particles (shown in Fig. 5). 
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Some of these fragments eventually merge with other more 
massive sinks, but by the end of our simulation we have 
5 sink particles throughout the disk, described in Table 1. 
The two largest sinks are 700 AU apart, have a mass ratio 
of ~ 1/3, and form 300 yr apart (sec Sec. 4.2.4 and Fig. 9). 
This is not unsimilar to the two fragments found in Turk 
et al. (2009), which had a separation of 800 AU, a mass ra- 
tio of ~ 3/5, and formed ~ 60 yr apart. The three smallest 
sinks in our calculation are all around 1 M© and formed after 
approximately 4000 yr. Compared to the two largest sinks, 
these smaller sinks are quite low-mass and have been present 
in the simulation for only a short time. The two largest sinks 
seem to be the only long-term sinks in the simulation. 

Analytical estimates concerning massive disk fragmen- 
tation by Kratter & Matzner (2006) showed interesting sim- 
ilarity to our results. They find that protostellar disks host- 
ing massive stars (O and B stars) are prone to fragmenation 
for disk radii larger than 150 AU, and they predict that 
many low-mass companions will form in the disks around 
these stars at initial separations of 100-200 AU (compare 
with our Table 1). Though their estimate is based on opti- 
cally thick disks undergoing viscous heating as well as stel- 
lar irradiation, the stabilizing effect of those heat sources is 
counteracted by rapid stellar accretion and high disk angu- 
lar momentum. If protostellar feedback were included in our 
work, this could reduce the number of fragments that form, 
since the initial protostar could heat the disk and inhibit 
further fragmentation, perhaps leaving a binary instead of 
a mutiple system. This effect is shown in Krumholz et al. 
(2007), where the evolution of high- mass prestellar cores in 
the present-day Universe is simulated. They find that pro- 
tostellar radiative feedback inhibits fragmentation and that 
the majority (<^ 90 percent) of accreted stellar mass goes 
onto one protostar, though a small number of lower-mass 
stars can form at sufficiently large distances from the main 
protostar or in disk fragments formed through Toomre in- 
stability. It will thus be important to include protostellar 
feedback in future simulations to fully ascertain the multi- 
plicity of Pop III star formation. 



4-2.4 Accretion rate 

How do the initial hydrostatic cores develop into fully 
formed Pop III stars? Fig. 9 shows the growth of the two 
most massive sink particles over time. We find that the main 
sink approximately grows as M oc t'^'^^ (red hne). By 5000 yr 
after sink formation, the main sink has grown to 43 Mq , only 
slightly higher than the corresponding result from Bromm & 
Loeb (2004). The jumps in mass arc instances where smaller 
sinks merged with the largest sink. One caveat to keep in 
mind is that we did not explicitly check that these smaller 
sinks were not centrifugally supported against infall towards 
the main sink before they were merged. However, we exam- 
ined the largest secondary sink (~ 4 M0) that merged with 
the main sink. This merger occured at 3700 yr after ini- 
tial sink formation. We checked the secondary sink's specific 
angular momentum at the last simulation output before it 
merged with the main sink. At this time it was ~ 100 AU 
away from the main sink, and its specific angular momen- 
tum with respect to the main sink was 2x10^^ cm^ s~^. This 
is less than the specific angular momentum required at this 



distance for centrifugal support against infall, 2.6 xlO^'^ cm^ 
. It is thus plausible that this was indeed a true merger. 

In Fig. 9, we also show the corresponding accretion 
rate onto the most massive sink particle over time (panel 
b). The accretion rate is highly variable and initially very 
large: M ~ 7cl/G, if we evaluate the sound speed at 700 K, 
and ~ 50cl/G at 200 K. These rates are similar to those 
predicted in analytic solutions (Larson 1969; Penston 1969; 
Shu 1977), and to the result found in Bromm & Loeb (2004). 
Note that the accretion rate does generally decline over time 
approximately as M oc t"" ''^, except for a few 'spikes' in ac- 
cretion when the main sink merges with other smaller sinks. 
For comparison the accretion rate onto the second largest 
sink (dotted line) is also shown, whose powerlaw fit goes as 
M oc t~'^'^^ (fit not shown) and is normalized to a value of 
about 1/5 that of the main sink. The lower rate is expected 
since the second largest sink accretes significantly cooler gas. 

As accretion onto the main sink continues, the overall 
average accretion rate will drop to a much lower value, closer 
to the average value found in Bromm & Loeb (2004). If 
accretion were to continue for the entire lifetime of the star, 
3 X 10® Myr, and if the accretion rate continued to decrease 
as approximately M oc t"" *^, the final Pop HI star would 
reach M« ^ 1000 Mq. However, our fit likely overestimates 
the accretion rate at late times, and the actual final mass 
of the star is likely to be significantly lower. Furthermore, 
feedback effects will likely slow accretion, or even halt it 
entirely, before the star dies. Thus, our extrapolation serves 
as a robust upper limit to the final Pop HI stellar mass. 

4.2.5 Thermodynamics of accretion flow 

After the initial sink particle has grown in mass, the sur- 
rounding gas divides into two phases - a hot and cold one. 
In Fig. 10, we illustrate this bifurcation with a temperature- 
density phase diagram at various stages of the simulation. 
Heating becomes significant once the initial sink grows be- 
yond 10 Mq . At this mass the gravitational force of the sink 
particle is strong enough to pull gas towards it with velocites 
sufficiently high to heat the gas to a maximum temperature 
of ~ 7,000 K (see discussion in Sec. 4.2.1). This is the tem- 
perature where the collisional excitation cooling of atomic 
hydrogen begins to dominate over the adiabatic and viscous 
heating. The increasing mass of the sink causes a pressure 
wave to propagate outward from the sink and heat particles 
at progressively larger radii and lower density, as is visible 
in Fig. 5 and Fig. 10. The heated region extends out to 
~ 2,000 AU at 5000 yr. The pressure wave thus propagates 
at a subsonic speed of ~ 2 km s^'^. 

The main sink is able to accrete a fraction of these 
heated particles (e.g. the yellow particle in Fig. 10), while 
it continues to accrete cold particles from the disk as well. 
We record the temperatures of the SPH particles accreted 
onto the sink, so that we can track the relative contribu- 
tions from the two phases. When the main sink has grown 
to ~ 30 Mq, accretion from the hot phase begins to occur, 
such that heated particles contribute slightly over 50% to the 
total rate towards the end of the simulation. By this time, 
the average temperature of all particles accreted onto the 
main sink is close to 3000 K, a value indicative of significant 
contributions from the hot phase. The second-most massive 
sink, on the other hand, accretes particles from only the 
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Figure 9. (a) Sink mass versus time since initial sink formation. The solid line shows the mass of the first sink particle over time in 
our calculation. The red line is the least-squares powerlaw fit to the sink mass over time, which goes as M oc ("-SS. The dash-dotted line 
shows the result obtained by Bromm & Loob (2004). The dotted line is the mass growth of the second largest sink, (b) Accretion rate 
versus time since initial sink formation. The solid line shows the accretion rate of the first sink particle throughout our calculation. The 
red line is the powerlaw fit to the mass accretion rate (M oc t"" *^). The black dash-dotted line shows the instantaneous accretion rate 
found by Bromm & Loeb (2004), while the red dotted line is their average accretion rate over the lifetime of a Pop III star, determined 
by extrapolating their rate to 3x10^ yr. The dotted line is the accretion rate of the second largest sink. Also shown are the Shu accretion 
rates for isothermal gas at 700 K (green short-dash line) and 200 K (blue long-dash line). 



cool phase as it has not yet reached a sufficiently high mass 
to capture heated particles. It has, however, just become 
massive enough to begin heating its surrounding particles. 
Fig. 11 reveals the location of the two diflerent gas phases 
in a complementary way by showing where cooling due to 
H2 (left panel) and atomic hydrogen line cooling dominate 
(right panel) . H2 cooling is important in the cold, dense gas 
of the disk-like structure, whereas atomic cooling occurs in 
the hot, nearly spherical region around the largest sink par- 
ticle. This region grows in size as the sink itself gains mass. 
By the end of the simulation, a similar heated region is also 
beginning to expand around the second-most massive sink. 

Since the heated region is created by the sink parti- 
cle's potential well, we can estimate its extent by evaluating 
the Bondi radius, which marks the distance out to where 
the gravitational energy of the sink will be greater than the 
thermal energy of the gas. If we use the final mass of the 
main sink, Msi^k — 40 Mq, and the maximum soundspeed 
of the gas particles, Cs ~ 7 km s~^ for temperatures of 7000 
K, we find 



Rb 



750AU 



(10) 



This is close to the size of the heated region seen in both 
Fig. 5 and Fig. 11. We should further note that, except for 
sink masses lower than ~ 2 M0, Rb, will be larger than race, 
which is held constant at 50 AU. When the main sink reaches 
larger masses and begins to accrete heated particles, Rb is in 
fact over an order of magnitude larger than race, indicating 
that our treatment of accretion is indeed conservative and 
that there is little un-physical accretion of heated particles. 



4.2.6 Feedback 

As noted earlier, the formation of a disk-like configuration 
around the protostar will have important implications for 
feedback effects. As demonstrated by McKee & Tan (2008) , 
the disk structure will mitigate the impact of protostellar 
feedback. Radiation will be able to escape along the polar 
directions, and although the innermost, optically thick, part 
of the accretion disk is not resolved in our simulation, it 
would shield a portion of the outer accretion flow from direct 
feedback. Future three-dimensional simulations that include 
radiative feedback will yield a better understanding of how 
and when disk accretion onto a primordial protostar may 
be terminated. We now wish to assess how important this 
neglected feedback would be up to the stage simulated here. 

In Fig. 12, we compare the accretion luminosity. Lace, 
with the Eddington luminosity due to electron scattering 
and due to H~ opacity. Determining Lace requires an esti- 
mate of the protostellar radius /?* . In the initial adiabatic ac- 
cretion phase, before Kelvin-Helmholtz contraction has com- 
menced, the photospheric opacity of the protostar is domi- 
nated by H~ bound- free absorption (i.e., 'Phase P of Omukai 
& Palla 2003). The extremely strong sensitivity of the H~ 
bound-free opacity to temperature (kjj- oc T^'^'^) locks the 
photospheric temperature to ^ 6000 K. We estimate that 
the protostar emits as a blackbody at this temperature and 
furthermore assume that the protostellar luminosity during 
this phase is dominated by Lace. This is justified as long as 
iace ^ tKK, where taec is the accretion timescale and Ikh the 
Kelvin-Helmholtz time. We then have 



L,7 



GM,M 



A-KRiiasBT"^ , 



(11) 



where L«/ is the protostellar luminosity during the adiabatic 
accretion phase, asB is the Stefan-Boltzmann constant, and 
T = 6000 K. This results in 
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Figure 10. Evolution of thermodynamic properties. We plot temperature vs. number density at various times during the simulation. 
The colored dots trace the position of four representative particles, illustrating how the hot phase is populated. Note the hot particle 
that is accreted by the main sink sim,2000yr after the sink's creation (yellow dot). The times are as follows: (a) 1000 yr. (b) 2000 yr. 
(c) 3000 yr. (d) 5000 yr. 




Figure 11. Dominant cooling processes in the center of the minihalo. Sinks are denoted as in Fig. 5. Left : H2 cooling rate within the 
central 5000 AU, shown in projection along the x-z plane, after ~ 5000 yr of accretion. Right : H line cooling rate in the same region and 
at the same time, again shown in projection along the x-z plane. 
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50Rp 



1/3 



M 



1/3 



(12) 



where ii*/ is the protostellar radius during the adiabatic 
accretion phase, and Mm ~ 4.4 x 10"'^ Mq yr~^, the fiducial 
accretion rate used by Stabler et al. (1986) and Omukai & 
Palla (2003). Our simple estimate of how varies with 
mass and accretion rate is quite similar to the results of these 
more detailed previous studies, and here serves to highlight 
the relevant physics. 

During the next phase, the protostellar radius begins 
to shrink due to Kelvin-Helmholtz contraction. The major 
source of opacity in the interior of the protostar is elec- 
tron scattering, which is independent of density and tem- 
perature. This along with hydrostatic equilibrium yields the 
mass-luminosity relation L» oc M'i . Using figure 4 of Omukai 
& Palla (2003), we normalize this relation as follows: 



lO^Lo 



M, 



10 M, 



(13) 



where L,// is the protostellar luminosity during the Kelvin- 
Hehnholtz contraction phase. If we further assume that tKH 
— ^acc, we arrive at 



GMj 



M 



(14) 



Using L« ~ -L*/7, we now find a relation for R^ii, the 
protostellar radius during the Kelvin-Helmholtz contraction 
phase: 



R,i 



140Rp 



(S-^( 



M, 



(15) 



In determining the evolution of Lace, we set 7?, — R,i at 
early times when the protostar is still experiencing adiabatic 
accretion. When M, ~ 10 Mq , after approximately 800 yr of 
accretion, the value of R,ii falls below that of R,i. At this 
point, we estimate that Kelvin-Helmholtz contraction has 
commenced and switch to setting R, — R,ii. The resulting 
transition in the growth of Lace can be seen in Fig. 12. The 
discontinuous slope in the accretion luminosity is an artefact 
of our approximate modelling of protostellar structure. 

Since the accretion luminosity does not exceed either 
Eddington luminosity, we conclude that for the duration of 
our simulation, radiation pressure is not likely to signifi- 
cantly change the early accretion rates we have found. At 
later times, when the protostar reaches higher masses, how- 
ever, and in particular when it reaches the ZAMS and hy- 
drogen burning turns on, such feedback effects will become 
more important, and can no longer be neglected. 



5 PARAMETER SENSITIVITY 

How commonly do Pop HI stars form in multiples? To what 
extent does this result depend on the initial cosmological pa- 
rameter values, particularly the value for erg? Is there a rela- 
tion between the spin of a minihalo and the probability that 
it will host a Pop HI multiple (see section 4.1)? To address 
this uncertainty, we make an initial estimate of the sensitiv- 
ity of our results to differing values of as by comparing the 
specific angular momentum profile of the central 1000 M0 of 
our simulation with that of Abel et al. (2002) and Yoshida 



10^ 



10' 



10- 



1000 2000 3000 4000 5000 
time [yr] 



Figure 12. Impact of feedback on the accretion process. Solid 
red line is the accretion luminosity Lace of the most massive sink 
versus time, calculated using i?« = R^j at early times and /?« = 
Rtll at later times, as described in the text. In determining Lace 
we use the power-law fit to the accretion rate (M oc t"'''*^). 
Dotted line is the prediction of Lace by Tan & Mckee (2004) given 
/d = 1 and K' = 1.33. Dashed line is the Eddington luminosity 
Lfidd due to electron scattering, calculated using the sink mass. 
The dash - doted line is Ljjdd due to H~ bound-free opacity, 
calculated at a temperature of 6000 K. Note that the kink in the 
curves for the accretion luminosity is an artefact of our idealized 
modelling of protostellar structure. 



et al. (2006) (Fig. 13). The angular momentum profile for 
our work is taken at the point immediately before the first 
sink forms. Abel et al. (2002) initialized their simulation us- 
ing h=0.5, f^A ~ 0, fie = 0.06 and erg = 0.7. Yoshida et al. 
(2006) used h=0.7, D.a = 0.7, SIb = 0.04 and as = 0.9. 
If the particular parameter values used in initializing our 
calculation made the central gas more prone to develop a 
flattened disk- like structure that would later fragment, this 
would first be apparent in an enhanced angular momentum 
profile. However, despite the varying initial conditions used 
for each of these calculations, the angular momentum pro- 
files are all quite similar in shape and amplitude, especially 
in the innermost region. The central 100 Mq that will later 
become the star-disk system in our calculation does not dif- 
fer in angular momentum from the other two calculations 
by more than a factor of 2 to 3. In fact, in this mass range 
the simulation of Yoshida et al. (2006) tends to have the 
highest angular momentum. The initial angular momentum 
distribution of the star-forming gas thus seems fairly robust 
within the range of cosmological parameter values sampled 
with these three simulations, and our overall result should 
be correspondingly robust. We further note that Clark et al. 
(2008) find that their non-cosmological simulation of a pri- 
mordial gas cloud also resulted in a profile very similar to 
the same two cosmological simulations, and upon extending 
their calculation through the sink particle method they also 
found extensive fragmentation. 

However, the generality of our result will be best de- 
termined through a more comprehensive set of simulations. 
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Figure 13. Comparison of specific angular momentum profiles 
for three different cosmological simulations: this work at the point 
immediately before initial sink formation (solid line), Yoshida et 
al. (2006, dotted Une), and Abel et al. (2002, dashed line). Despite 
the differing initial conditions for each simulation, the initial an- 
gular momentum distributions are quite similar, particularly at 
the smallest enclosed mass scales from which most of the mass of 
the sinks are accreted. 

aiming at studying fragmentation over a range of minihalo 
spins. We will report on the results in a future work. Fi- 
nally, wc expect that the effects from protostellar feedback, 
neglected here, will likely dominate over those arising from 
variations in the initial conditions. Again, this needs to be 
addressed with improved simulations. 



6 SUMMARY AND DISCUSSION 

We have evolved a cosmological simulation until a primor- 
dial minihalo has formed at z ~ 20, and follow the sub- 
sequent collapse of the gas in its center up to densities of 
n = lO^^cm"^. We find that the initial collapse, leading 
to the formation of a small hydrostatic core, is very similar 
to previous high-resolution work. Subsequently, a disk-like 
structure grows around the central core, eventually compris- 
ing a mass of ~ 40 . Wc use the sink particle method to 
study the mass flow onto the central protostar for 5000 yr af- 
ter its formation. Similar to the fragmentation found in Turk 
et al. (2009), wc find that a second protostar forms soon after 
the first one. Upon further evolution, the disk exhibits spi- 
ral structure and undergoes further fragmentation, resulting 
in several more sink particles by the end of the simulation. 
Most of the mass of the smaller, secondary sinks is accreted 
through the disk, which consists mostly of cold gas. Once the 
main sink grows beyond 10 M0, its strong gravity is able to 
heat the surrounding particles to several thousand Kelvin, 
such that a two-phase medium emerges. The main sink is 
able to accrete these heated particles as well. After 5000 yr 
of accretion, the most massive sink has reached ^^40 Mq, at 
which point we terminate our sirrmlation. Beyond this time, 
protostellar feedback effects, which are not included here, 
could no longer be neglected. The Pop III star, however, is 



expected to accumulate more mass. Even in the no-feedback 
case, however, the finite time available for accretion limits 
the Pop III mass to ^ 1000 M©, where it is evident that in 
realistic settings, masses will not reach such extreme values. 
The accretion rate and the rmmber of stars that form will 
be affected by protostellar feedback, with the likely outcome 
that both individual masses and the number of stars will be 
reduced (see Krumholz et al. 2007). 

The recent radiation-hydrodynamical simulation by 
Krumholz ct al. (2009) , which studies the collapse of a mas- 
sive pre-stellar core in a present-day GMC, shows interesting 
similarity to our investigation of the primordial case. They 
followed the collapse of a 100 M0 gas cloud, the formation 
of several protostellar seeds and their subsequent growth 
by accretion. Like our simulation, their system also evolved 
into a disk with a radius of order 1000 AU which fragmented 
due to gravitational instability. Similar to the sink mergers 
seen in our calculation, in their simulation the secondary 
stars that formed in the disk collided with the most massive 
star. By the end of their simulation, which was halted af- 
ter ^ 60,000 yr, the system had evolved into a binary with 
stars of mass ~ 30 and 40 Mq , masses fairly close to those 
of the two largest sinks in our simulation. Their calculation, 
however, included dust opacity as well as radiation pressure 
from the stars, and their initial conditions were typical of 
present-day stax formation. Their accretion and fragmenta- 
tion timescales were also much longer than those found in 
our calculation, but the overall qualitative results were alike. 

Thus, the formation of binary and multiple systems 
within disk structures is likely not limited to modern-day 
star formation and can be found even in the primordial case, 
arising from cosmological initial conditions. Binary and rrml- 
tiple systems may be much more common in the Pop III case 
than previously thought. As similarly argued in Clark et al. 

(2008) , this is because most cosmological simulations do not 
follow the gas evolution for many years after the first proto- 
star has formed, and fragmentation may not typically occur 
until after the simulations are stopped. For instance, the 
second protostar in our calculation did not form until 300 yr 
after the first, while the calculation of Turk et al. (2009) was 
stopped only 200 yr after the first protostar formed. 

If Pop III binaries and multiples are actually com- 
mon, then this has important implications for the final fate 
and observational signature of Pop III.l stars, those zcro- 
metallicity stars that have not yet been influenced by any 
previous stellar generation (e.g. Bromm et al. 2009). The 
presence of multiple accretors in our simulation did not seem 
to affect the accretion rate onto the largest sink as compared 
to previous work (Bromm & Locb 2004). However, extend- 
ing our simulation to later times might show that the most 
massive stax would ultimately have to compete with its sec- 
ondary companion for mass. As pointed out in Turk et al. 

(2009) , if this in the end prevents the largest star from reach- 
ing very high masses (^ 100 M0), this would also affect its 
resulting chemical signature. If multiple star formation were 
common in Pop III, then the occurrence of pair-instability 
supernovae (PISNe) might be significantly reduced, since 
stars that explode as PISNe must be in the higher mass 
range (140 Mq < M, < 260 Mq, Heger & Woosley 2002). 
Extremely metal-poor Galactic halo stars are believed to 
preserve the nucleosynthetic pattern of the first SNe, so 
this may help to explain the lack of PISNe chemical sig- 
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natures found in halo stars (e.g. Christlieb at al. 2002; Beers 
& Christlieb 2005; Frebel et al. 2005; Tumlinson 2006; Karls- 
son et al. 2008). Turk et al. (2009) also note that if multi- 
ple high-mass Pop III stars can form within a minihalo at 
similar times, the ionizing flux on neighboring minihaloes 
could be increased, possibly leading to more efhcient molec- 
ular cooling and ultimately lower-mass (Pop III. 2) stars in 
the neighboring haloes (e.g. O'Shea et al. 2005). This again 
would influence the chemical signature of extremely metal- 
poor Galactic halo stars. 

It should be mentioned that tight binaries on scales 
smaller than the resolution of our simulation could also form. 
Though a higher resolution study would bo necessary to con- 
firm this, the formation of a fragmenting disk of steadily in- 
creasing angular momentum in our cosmological simulation 
makes this a reasonable possibility. This could also limit the 
mass of Pop III stars, yielding similar effects on their re- 
sulting chemical signature as discussed above. The possible 
binary nature of Pop III stars, provided that they would not 
grow to become very massive, may also yield a direct pri- 
mordial formation pathway to carbon-enhanced metal-poor 
(CEMP) stars, those metal-poor halo stars with unusually 
high carbon abundance (e.g. Frebel ct al. 2007; Tumlinson 
2007). If one of the members of a Pop III binary was an 
intermediate-mass star (1 M© ^ M, <, 8 M©), then the 
companion might undergo carbon enhancement through bi- 
nary mass-transfer (e.g. Suda et al. 2004). This would occur 
when the intcrrncdiatc-mass star enters an asymptotic giant 
branch (AGB) phase during which its C-rich ejecta can be 
captured by the companion. Furthermore, another type of 
Pop III signature would arise if gamma-ray bursts (GRBs) 
were able to form from very tight Pop III binaries. These 
may be detectable by Swift (see Bromm & Loeb 2006; Bel- 
czynski et al. 2007). 

On the other hand, if both members of a typical Pop III 
binary were in the mass range to collapse directly into a 
black hole (M« ^ 260 Mq), then they could become mas- 
sive black hole (MBH) binaries. Such a population of MBH 
binaries may emit a gravitational wave (GW) signature de- 
tectable by the future space-born Laser Interferometer Space 
Antenna (LISA). Previous studios have already explored the 
possibility of GW detection of MBH binaries formed through 
halo mergers during early structure formation (e.g. Wyithe 
& Loeb 2003; Sesana et al. 2005). Sesana et al. (2005) find 
that high redshift {z ^ 10) MBH binaries with masses ^ 10^ 
Mq will be detectable by LISA. MBH binaries formed from 
Pop III progenitors thus approach the range of detectability 
by LISA, especially if the BHs continue to grow in mass 
through accretion. Furthermore, Belczynski et al. (2004) 
predict that advanced LIGO could also detect the GWs 
from Pop HI MBH binary mergers. Madau & Rees (2001) 
suggested a different type of Pop III MBH signal when they 
discussed how an accreting primordial MBH could tidally 
capture a star and become a detectable off-nuclear ultralu- 
minous X-ray source. They conclude, however, that such a 
non- primordial origin to Pop HI binaries is extremely inef- 
ficient. 

A modified picture of Pop HI star formation is thus 

emerging in which disk formation, fragmentation, and the 
resulting binary and mutiplc systems will play an important 
role in determining the mass of Pop HI stars and their in- 
fluence on later star and galaxy formation. Further studies. 



particularly ones which include protostellar feedback, will 
continue to add to this developing picture. 
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